Santa Barbara County Vineyards Analysis 2023

Authors:
Gregory Halverson and Claire Villanueva-Weeks
Jet Propulsion Laboratory, California Institute of Technology

This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, and was sponsored by ECOSTRESS and the National Aeronautics and Space Administration (80NM0018D0004).

© 2023. All rights reserved.

Summary and Background

Santa Barbara County is a well-known for its viticulture. Positioned near the ocean, the unique topography of Santa Barbara allows for unique microclimates that yield diverse crops and wines. In this notebook, we aim to visualize and analyze Vineyard health using Earth observation data from two different satellite products: the EMIT L2A Reflectance product stored in NETCDF4 format and an ECOSTRESS L2 Land Surface Temperature product stored in GEOTIFF format. Our goal is to visualize Normalized Difference Vegetation Index (NDVI), Equivalent Water Thickness (EWT), Evapotranspiration (ET), and Surface Temperature (ST) interactively in holoviews maps and generate linear regression plots of point source data.

Figure 1. Vineyard.

Requirements

• NASA Earthdata Account
  
• No Python setup requirements if connected to the workshop cloud instance!
  
• Local Only Set up Python Environment - See setup_instructions.md in the /setup/ folder to set up a local compatible Python environment
  
• Downloaded necessary files. This is done at the end of the 01_Finding_Concurrent_Data notebook.

Tutorial references

01_Finding_Concurrent_Data.ipynb - Finding and downloading concurrent EMIT and ECOSTRESS data with earthaccess

03_EMIT_CWC_from_Reflectance.ipynb - Calculating EWT of ROI

03_EMIT_Inreractive_Points.ipynb - Generating a map with clickable points

Required datasets

• EMIT Reflectance
  • EMIT_L2A_RFL_001_20230401T203803_2309114_003.nc
• EMIT Canopy Water Content
  • EMIT_L2A_RFL_001_20230401T203751_2309114_002_roi_bbox_cwc_merged.tif
• ECOSTRESS Land Surface Temperature
  • ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST.tif
• ECOSTRESS Evapotranspiration
  • ECOv002_L3T_JET_26860_001_10SGD_20230401T203732_0700_01_ETdaily.tif
• SB_ROI_ag.geojson: This file includes polygons of delineated agricultural field boundaries. It has been clipped to our ROI and fields were edited for ease of usability.

Tutorial Outline

5.1 Set up
5.2 Locating Data Sources and Loading Data
5.3 Color Hex-Codes for Raster Plots
5.4 Subsetting the data over our ROI
5.5 Interactive Maps
5.6 Scatter Plots

References

1. Jackson, M. (2018). Santa Barbara County (Images of America). Arcadia Publishing.
2. CalAgData

5.1 Set up

These are some built-in Python functions we need for this notebook, including functions for handling filenames and dates.

We’re using the rioxarray package for loading raster data from a GeoTIFF file, and we’re importing it as rxr. We’re using the numpy library to handle arrays, and we’re importing it as np. We’re using the rasterio package to subset the data.

We’re using the geopandas library to load vector data from GeoJSON files, and we’re importing it as gpd. We’re using the shapely library to handle vector data and the pyproj library to handle projections.

Import the emit_tools module and call use from emit_tools import emit_xarray help(emit_xarray) the help function to see how it can be used.

Note: This function currently works with L1B Radiance and L2A Reflectance Data.

# Some cells may generate warnings that we can ignore. Comment below lines to see.
import warnings
warnings.filterwarnings('ignore')

import os
from os.path import join, expanduser, splitext, basename
import numpy as np
import pandas as pd
import geopandas as gpd
from osgeo import gdal
import rasterio as rio
import rioxarray as rxr
import holoviews as hv
import hvplot.xarray
import hvplot.pandas
import geoviews as gv
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

from holoviews.streams import SingleTap
from scipy.stats import linregress

Importing Custom EMIT Tools to handle the EMIT reflectance data

from modules.emit_tools import emit_xarray, ortho_xr

We can define these constants to prescribe the dimensions of our figures. Feel free to adjust these to fit your display. We’re setting the alpha to make the raster semi-transparent on top of the basemap.

FIG_WIDTH_PX = 1080
FIG_HEIGHT_PX = 720
FIG_WIDTH_IN = 16
FIG_HEIGHT_IN = 9
FIG_ALPHA = 0.7
BASEMAP = "EsriImagery"
SEABORN_STYLE = "whitegrid"
FILL_COLOR = "none"
LINE_COLOR = "red"
LINE_WIDTH = 3
WEBMAP_PROJECTION = "EPSG:4326"
sns.set_style(SEABORN_STYLE)

5.2 Locating Data Sources and Loading Data

Here we are defining location of EMIT and ECOSTRESS raster files.

data_dir = "../data/" 
data_dir

'../data/'

Here we are loading in our vector data source. We can image the delineated agricultural field boundaries and take a look at the fields.

ag_file = join(data_dir,"SB_ROI_ag.geojson")
# cut out the agroi section 
ag_latlon = gpd.read_file(ag_file)

ag_fig = ag_latlon.to_crs(WEBMAP_PROJECTION).hvplot.polygons(
    tiles=BASEMAP,
    line_color = LINE_COLOR,
    line_width = LINE_WIDTH,
    fill_color = FILL_COLOR,
    crs='EPSG:4326'
) 

ag_fig

Unable to display output for mime type(s): 

ag_latlon

	OBJECTID	permittee	site_name	size	size_units	crop_list	calc_acres	mtrs	geometry
0	3654110	CENT ANNI-COASTAL VINEYARD CARE	CENT ANNI	8.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	1.463097	S07N34W11	MULTIPOLYGON (((-8064.151 -373076.101, -8058.8...
1	3671186	CENT ANNI-COASTAL VINEYARD CARE	CENT ANNI	8.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	1.232124	S07N34W11	MULTIPOLYGON (((-8220.677 -373296.021, -8205.1...
2	3671187	CENT ANNI-COASTAL VINEYARD CARE	CENT ANNI	8.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	2.782784	S07N34W11	MULTIPOLYGON (((-8185.460 -372957.510, -8146.0...
3	4829810	CHUMASH VINEYARDS LLC- COASTAL VINEYARD	CAMP 4 -	280.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	277.792034	S06N30W05	MULTIPOLYGON (((-5734.564 -375663.536, -5340.6...
4	3681606	GREAT OAKS RANCH-FINKLE	GREAT OAKS RANCH VINEYARD - ALL PLOTS	20.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	15.713672	S07N30W30	MULTIPOLYGON (((-8104.976 -374002.479, -8098.9...
...	...	...	...	...	...	...	...	...	...
463	8928511	BROTHERS BEST PRODUCE	ARBELAITZ RANCH	12.0	ACRES	OAT, Rotational Crops	11.131229	S07N31W36	MULTIPOLYGON (((-8961.836 -375463.076, -8960.0...
464	8928512	BROTHERS BEST PRODUCE	SITE 06: MORO RANCH TRI-CLOR EC	15.0	ACRES	SQUASH	14.690810	S07N30W32	MULTIPOLYGON (((-4728.422 -374455.107, -4732.8...
465	8928514	BROTHERS BEST PRODUCE	SHEPPARD RANCH	108.0	ACRES	OAT, Rotational Crops	110.279936	S06N30W05	MULTIPOLYGON (((-6938.727 -376457.113, -5628.6...
466	8945655	CAMARILLO FARMING	SITE 01 BASELINE AVE RANCH	8.0	ACRES	PEPPER FRUITNG, TOMATILLO, CABBAGE	7.910668	S07N31W36	MULTIPOLYGON (((-8567.652 -375377.567, -8567.0...
467	8947580	TANNER VINEYARD	TANNER RANCH	1.0	ACRES	GRAPE, WINE	1.006459	S06N31W01	MULTIPOLYGON (((-8856.269 -376329.948, -8907.5...

468 rows × 9 columns

Loading in EMIT reflectance data

To open up the EMIT .nc file we will use the netCDF4, xarray and emit_tools libraries.

EMIT_fp = join(data_dir, "EMIT_L2A_RFL_001_20230401T203803_2309114_003.nc") 
EMIT_fp

'../data/EMIT_L2A_RFL_001_20230401T203803_2309114_003.nc'

refl_ds = emit_xarray(EMIT_fp, ortho=True)
refl_ds

<xarray.Dataset>
Dimensions:           (latitude: 1895, longitude: 2368, wavelengths: 285)
Coordinates:
  * wavelengths       (wavelengths) float32 381.0 388.4 ... 2.486e+03 2.493e+03
    fwhm              (wavelengths) float32 ...
    good_wavelengths  (wavelengths) float32 ...
  * latitude          (latitude) float64 35.31 35.3 35.3 ... 34.28 34.28 34.28
  * longitude         (longitude) float64 -120.3 -120.3 -120.3 ... -119.1 -119.1
    elev              (latitude, longitude) float32 nan nan nan ... nan nan nan
    spatial_ref       int64 0
Data variables:
    reflectance       (latitude, longitude, wavelengths) float32 nan nan ... nan
Attributes: (12/40)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-1.20348784e+02  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Estimated Surface Reflectance...
    granule_id:                        EMIT_L2A_RFL_001_20230401T203803_23091...
    Orthorectified:                    True

xarray.Dataset

• Dimensions:
  • latitude: 1895
  • longitude: 2368
  • wavelengths: 285
• Coordinates: (7)
  • wavelengths
    (wavelengths)
    float32
    381.0 388.4 ... 2.486e+03 2.493e+03

    long_name :

    Wavelength Centers

    units :

    nm

    array([ 381.00558,  388.4092 ,  395.81583, ..., 2478.153  , 2485.5386 ,
           2492.9238 ], dtype=float32)

  • fwhm
    (wavelengths)
    float32
    ...

    long_name :

    Full Width at Half Max

    units :

    nm

    [285 values with dtype=float32]

  • good_wavelengths
    (wavelengths)
    float32
    ...

    long_name :

    Wavelengths where reflectance is useable: 1 = good data, 0 = bad data

    units :

    unitless

    [285 values with dtype=float32]

  • latitude
    (latitude)
    float64
    35.31 35.3 35.3 ... 34.28 34.28

    long_name :

    Latitude (WGS-84)

    units :

    degrees north

    array([35.305249, 35.304707, 35.304165, ..., 34.279345, 34.278803, 34.278261])

  • longitude
    (longitude)
    float64
    -120.3 -120.3 ... -119.1 -119.1

    long_name :

    Longitude (WGS-84)

    units :

    degrees east

    array([-120.348513, -120.347971, -120.347428, ..., -119.066133, -119.065591,
           -119.065049])

  • elev
    (latitude, longitude)
    float32
    nan nan nan nan ... nan nan nan nan

    long_name :

    Surface Elevation

    units :

    m

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]], dtype=float32)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Data variables: (1)
  • reflectance
    (latitude, longitude, wavelengths)
    float32
    nan nan nan nan ... nan nan nan nan

    long_name :

    Surface Reflectance

    units :

    unitless

    array([[[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
    ...
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]],
    
           [[nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            ...,
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan],
            [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)

• Indexes: (3)
  • wavelengths
    PandasIndex

    PandasIndex(Float64Index([ 381.0055847167969,  388.4092102050781,  395.8158264160156,
                  403.22540283203125, 410.63800048828125,  418.0535888671875,
                   425.4721374511719,  432.8927001953125,  440.3172607421875,
                   447.7427978515625,
                  ...
                   2426.440185546875,  2433.830322265625,   2441.21826171875,
                     2448.6064453125,  2455.994384765625,  2463.381591796875,
                   2470.767822265625,  2478.153076171875,   2485.53857421875,
                      2492.923828125],
                 dtype='float64', name='wavelengths', length=285))

  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 35.30524898855477,  35.30470675603452,  35.30416452351426,
                   35.30362229099401, 35.303080058473746,  35.30253782595349,
                   35.30199559343323,  35.30145336091298,  35.30091112839272,
                   35.30036889587247,
                  ...
                   34.28314068787152,  34.28259845535126,  34.28205622283101,
                   34.28151399031075, 34.280971757790496, 34.280429525270236,
                   34.27988729274998,  34.27934506022973,  34.27880282770947,
                   34.27826059518922],
                 dtype='float64', name='latitude', length=1895))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([-120.34851295716088, -120.34797072464062, -120.34742849212037,
                  -120.34688625960011, -120.34634402707985, -120.34580179455959,
                  -120.34525956203935, -120.34471732951909, -120.34417509699882,
                  -120.34363286447856,
                  ...
                  -119.06992867439637,  -119.0693864418761, -119.06884420935585,
                   -119.0683019768356, -119.06775974431534, -119.06721751179508,
                  -119.06667527927483, -119.06613304675457, -119.06559081423431,
                  -119.06504858171405],
                 dtype='float64', name='longitude', length=2368))

• Attributes: (40)

  ncei_template_version :

  NCEI_NetCDF_Swath_Template_v2.0

  summary :

  The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing (warming or cooling) of the atmosphere.\n\nThis file contains L2A estimated surface reflectances and geolocation data. Reflectance estimates are created using an Optimal Estimation technique - see ATBD for details. Reflectance values are reported as fractions (relative to 1). Geolocation data (latitude, longitude, height) and a lookup table to project the data are also included.

  keywords :

  Imaging Spectroscopy, minerals, EMIT, dust, radiative forcing

  Conventions :

  CF-1.63

  sensor :

  EMIT (Earth Surface Mineral Dust Source Investigation)

  instrument :

  EMIT

  platform :

  ISS

  institution :

  NASA Jet Propulsion Laboratory/California Institute of Technology

  license :

  https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/

  naming_authority :

  LPDAAC

  date_created :

  2023-04-06T02:26:51Z

  keywords_vocabulary :

  NASA Global Change Master Directory (GCMD) Science Keywords

  stdname_vocabulary :

  NetCDF Climate and Forecast (CF) Metadata Convention

  creator_name :

  Jet Propulsion Laboratory/California Institute of Technology

  creator_url :

  https://earth.jpl.nasa.gov/emit/

  project :

  Earth Surface Mineral Dust Source Investigation

  project_url :

  https://emit.jpl.nasa.gov/

  publisher_name :

  NASA LPDAAC

  publisher_url :

  https://lpdaac.usgs.gov

  publisher_email :

  lpdaac@usgs.gov

  identifier_product_doi_authority :

  https://doi.org

  flight_line :

  emit20230401t203803_o09114_s000

  time_coverage_start :

  2023-04-01T20:38:03+0000

  time_coverage_end :

  2023-04-01T20:38:15+0000

  software_build_version :

  010610

  software_delivery_version :

  010610

  product_version :

  V001

  history :

  PGE Run Command: {python /beegfs/store/emit/ops/repos/emit-sds-l2a/spectrum_quality.py /tmp/emit/ops/emit20230401t203803_emit.L2AReflectance_20230405t233102/output/emit20230401t203803_rfl /tmp/emit/ops/emit20230401t203803_emit.L2AReflectance_20230405t233102/output/emit20230401t203803_rfl_quality.txt}, PGE Input Files: {radiance_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203803/l1b/emit20230401t203803_o09114_s000_l1b_rdn_b0106_v01.img, pixel_locations_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203803/l1b/emit20230401t203803_o09114_s000_l1b_loc_b0106_v01.img, observation_parameters_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203803/l1b/emit20230401t203803_o09114_s000_l1b_obs_b0106_v01.img, surface_model_config=/beegfs/store/emit/ops/repos/emit-sds-l2a/surface/surface_20221020.json}

  crosstrack_orientation :

  as seen on ground

  easternmost_longitude :

  -119.06477746545393

  northernmost_latitude :

  35.3055201048149

  westernmost_longitude :

  -120.348784073421

  southernmost_latitude :

  34.277989478929086

  spatialResolution :

  0.000542232520256367

  spatial_ref :

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

  geotransform :

  [-1.20348784e+02 5.42232520e-04 -0.00000000e+00 3.53055201e+01 -0.00000000e+00 -5.42232520e-04]

  day_night_flag :

  Day

  title :

  EMIT L2A Estimated Surface Reflectance 60 m V001

  granule_id :

  EMIT_L2A_RFL_001_20230401T203803_2309114_003

  Orthorectified :

  True

Here some functions are set up to

5.3 Color Hex-Codes for Raster Plots

Here we have set up convenience functions for interpolating Color Hex-Codes for Raster Plots

def interpolate_hex(hex1, hex2, ratio):
    rgb1 = [int(hex1[i:i+2], 16) for i in (1, 3, 5)]
    rgb2 = [int(hex2[i:i+2], 16) for i in (1, 3, 5)]
    rgb = [int(rgb1[i] + (rgb2[i] - rgb1[i]) * ratio) for i in range(3)]
    
    return '#{:02x}{:02x}{:02x}'.format(*rgb)

def create_gradient(colors, steps):
    gradient = []
    
    for i in range(len(colors) - 1):
        for j in range(steps):
            ratio = j / float(steps)
            gradient.append(interpolate_hex(colors[i], colors[i+1], ratio))
    
    gradient.append(colors[-1])
    
    return gradient

def plot_cmap(cmap):
    gradient = np.linspace(0, 1, 256)  # Gradient from 0 to 1
    gradient = np.vstack((gradient, gradient))  # Make 2D image

    # Display the colormap
    plt.figure(figsize=(6, 2))
    plt.imshow(gradient, aspect='auto', cmap=cmap)
    plt.axis('off')
    plt.show()

Near-Infrared and Red Bands in EMIT Hyperspectral Cube

These are the nearest to 800 nm and 675 nm wavelengths, they can be used to calculate the NDVI using a ratio of the difference between the wavelengths to the sum of the wavelengths. NDVI is a metric by which we can estimate vegetation greenness.

The hvplot package extends xarray to allow us to plot maps. We’re reprojecting the raster geographic projection EPSG 4326 to overlay on the basemap with a latitude and longitude graticule. We will be using hvplot a few more times to look at the data we are using.

Defining Color-Map for Near Infrared

NIR_colors = ["#000000", "#FF0000"]
NIR_gradient = create_gradient(NIR_colors, 100)
NIR_cmap = LinearSegmentedColormap.from_list(name="NIR", colors=NIR_colors)
plot_cmap(NIR_cmap)

Extracting 800 nm Near Infrared from EMIT

NIR = refl_ds.sel(wavelengths=800, method="nearest")
NIR.rio.to_raster(join(data_dir,"NIR_800.tif"))
NIR = rxr.open_rasterio(join(data_dir,"NIR_800.tif")).squeeze("band", drop=True)
NIR_vmin, NIR_vmax = np.nanquantile(np.array(NIR), [0.02, 0.98])
NIR_title = f"{splitext(basename(EMIT_fp))[0]} ~800 nm"

NIR.rio.reproject(WEBMAP_PROJECTION).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap=NIR_gradient, 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(NIR_vmin, NIR_vmax),
    title=NIR_title
)

Unable to display output for mime type(s): 

Defining Color-Map for Red

red_colors = ["#000000", "#FF0000"]
red_gradient = create_gradient(red_colors, 100)
red_cmap = LinearSegmentedColormap.from_list(name="red", colors=red_colors)
plot_cmap(red_cmap)

Extracting 675 nm Red from EMIT

red = refl_ds.sel(wavelengths=675, method="nearest")
red.rio.to_raster(join(data_dir,"red_675.tif"))
red = rxr.open_rasterio(join(data_dir,"red_675.tif")).squeeze("band", drop=True)
red_vmin, red_vmax = np.nanquantile(np.array(red), [0.02, 0.98])
red_title = f"{splitext(basename(EMIT_fp))[0]} ~675 nm"

red.rio.reproject(WEBMAP_PROJECTION).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap=red_gradient, 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(red_vmin, red_vmax),
    title=red_title
)

Unable to display output for mime type(s): 

Defining Color-Map for Vegetation Index

NDVI_colors=[
    "#0000ff",
    "#000000",
    "#745d1a",
    "#e1dea2",
    "#45ff01",
    "#325e32"
]

NDVI_gradient = create_gradient(NDVI_colors, 100)
NDVI_cmap = LinearSegmentedColormap.from_list(name="NDVI", colors=NDVI_colors)
plot_cmap(NDVI_cmap)

Calculating Vegetation Index

NDVI = (NIR - red)/(NIR + red)
NDVI.rio.to_raster(join(data_dir,"NDVI.tif"))
NDVI = rxr.open_rasterio(join(data_dir,"NDVI.tif")).squeeze("band", drop=True)
NDVI_vmin, NDVI_vmax = np.nanquantile(np.array(NDVI), [0.02, 0.98])
NDVI_title = "Normalized Difference Vegetation Index"

NDVI_scene_map = NDVI.rio.reproject(WEBMAP_PROJECTION).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap=NDVI_gradient, 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(NDVI_vmin, NDVI_vmax),
    title=NDVI_title
)

NDVI_scene_map

Unable to display output for mime type(s): 

Loading in and visualizing our EWT data

We can now open up the EWT data that you calculated earlier with the LP DAAC EWT. The dataset is now in .tif format which means we can just open this one up

EWT_filename = join(data_dir,'EMIT_L2A_RFL_001_20230401T203751_2309114_002_roi_bbox_cwc_merged.tif')
EWT_filename

'../data/EMIT_L2A_RFL_001_20230401T203751_2309114_002_roi_bbox_cwc_merged.tif'

Defining Color-Map for EWT

EWT_colors = ["#FFFFFF", "#0000FF"]
EWT_gradient = create_gradient(EWT_colors, 100)
EWT_cmap = LinearSegmentedColormap.from_list(name="EWT", colors=EWT_colors)
plot_cmap(EWT_cmap)

Mapping EWT with color ramp

EWT = rxr.open_rasterio(EWT_filename).squeeze("band", drop=True)
EWT_vmin, EWT_vmax = np.nanquantile(np.array(EWT), [0.02, 0.98])
EWT_title = "Equivalent Water Thickness "

EWT.rio.reproject(WEBMAP_PROJECTION,nodata=np.nan).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap=EWT_gradient, 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(EWT_vmin, EWT_vmax),
    title=EWT_title
)

Unable to display output for mime type(s): 

Loading an ECOSTRESS LST granule

ST_filename = join(data_dir, "ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST.tif") 
ST_filename

'../data/ECOv002_L2T_LSTE_26860_001_10SGD_20230401T203733_0710_01_LST.tif'

The temperatures in the L2T_LSTE product are given in Kelvin. To convert them to Celsius, we subtract 273.15.

ST_K = rxr.open_rasterio(ST_filename).squeeze('band', drop=True)
ST_C = ST_K - 273.15

ST_C.rio.to_raster(join(data_dir,"ST.tif"))

ST_vmin, ST_vmax = np.nanquantile(np.array(ST_C), [0.02, 0.98])
ST_title = "ECOSTRESS Surface Temperature (Celsius)"

ST_C.rio.reproject(WEBMAP_PROJECTION,nodata=np.nan).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap="jet", 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(ST_vmin, ST_vmax),
    title=ST_title
)

Unable to display output for mime type(s): 

Defining Color Ramps

Color ramps, are used to map valued over a range of colors based on the data and goals of visualization. We’ve traditionally used the “jet” rainbow color-scheme to represent temperature, but this color-scheme is not accessible for all vision types.

Red-green colorblindness or deuteranomaly affects approximately 10% of people, so we are working on making a color ramp with that in mind. You can use hex codes to program in your own colors!

ST_colors = [
    "#054fb9",
    "#0073e6",
    "#8babf1",
    "#cccccc",
    "#e1ad01",
    "#f57600",
    "#c44601"
]

ST_gradient = create_gradient(ST_colors, 100)
ST_cmap = LinearSegmentedColormap.from_list(name="ST", colors=ST_colors)
plot_cmap(ST_cmap)

ST_vmin, ST_vmax = np.nanquantile(np.array(ST_C), [0.02, 0.98])
ST_title = "ECOSTRESS Surface Temperature (Celsius)"

ST_C.rio.reproject(WEBMAP_PROJECTION).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap=ST_gradient, 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(ST_vmin, ST_vmax),
    title=ST_title
)

Unable to display output for mime type(s): 

Loading ECOSTRESS Evapotranspiration granule

ET_filename = join(data_dir, "ECOv002_L3T_JET_26860_001_10SGD_20230401T203732_0700_01_ETdaily.tif")
ET_filename

'../data/ECOv002_L3T_JET_26860_001_10SGD_20230401T203732_0700_01_ETdaily.tif'

We have an established color-map for representing ECOSTRESS evapotranspiration that we believe to be inclusive.

ET_colors = [
    "#f6e8c3",
    "#d8b365",
    "#99974a",
    "#53792d",
    "#6bdfd2",
    "#1839c5"
]

ET_gradient = create_gradient(ET_colors, 100)
ET_cmap = LinearSegmentedColormap.from_list(name="ET", colors=ET_colors)
plot_cmap(ET_cmap)

The daily evapotranspiration product from ECOSTRESS is given in millimeters per day.

ET = rxr.open_rasterio(ET_filename).squeeze('band', drop=True)
ET_vmin, ET_vmax = np.nanquantile(np.array(ET), [0.02, 0.98])
ET_title = "ECOSTRESS Evapotranspiration (mm / day)"

ET.rio.reproject(WEBMAP_PROJECTION).hvplot.image(
    geo=True,
    rasterize=True, 
    cmap=ET_gradient, 
    alpha=FIG_ALPHA, 
    tiles=BASEMAP,
    clim=(ET_vmin, ET_vmax),
    title=ET_title
)

Unable to display output for mime type(s): 

5.4. Subsetting the data over our ROI

To clip the raster image to the extent of the vector dataset, we want to subset the raster to the bounds of the vector dataset. This dataset is included here in GeoJSON format, which we’ll load in as a geodataframe using the geopandas package.

The GeoJSON is provided by CalAg and includes information on agricultural field boundaries in the Santa Barbara County. We can use this code to select the largest polygons that are classified as vineyards.

Using the rioxarray clip function, we can subset the data using a geojson. to visualize the area that are vineyards are located, I have selected a boundary polygon. Overlaid are the outlines of some of the larger vineyards in the scene.

target_text = 'WINE'

filtered_gdf = ag_latlon[ag_latlon['crop_list'].str.contains(target_text, case=False, na=False)]

Q1 = filtered_gdf['size'].quantile(0.25)
Q3 = filtered_gdf['size'].quantile(0.75)
IQR = Q3 - Q1

# Define the threshold as Q3 + 1.5 * IQR (adjust multiplier as needed)
threshold_size = Q3 + 1.5 * IQR

# Filter polygons based on size greater than the threshold
vineyard_polygons = filtered_gdf[filtered_gdf['size'] > threshold_size]
vineyard_polygons

	OBJECTID	permittee	site_name	size	size_units	crop_list	calc_acres	mtrs	geometry
3	4829810	CHUMASH VINEYARDS LLC- COASTAL VINEYARD	CAMP 4 -	280.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	277.792034	S06N30W05	MULTIPOLYGON (((-5734.564 -375663.536, -5340.6...
53	3816318	PARKER STATION VINEYARDS-COASTAL VINEYAR	CAMP 4 - ALL PLOTS	260.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	25.585707	S06N30W05	MULTIPOLYGON (((-4888.679 -376414.168, -4885.6...
54	3791776	PARKER STATION VINEYARDS-COASTAL VINEYAR	CAMP 4 - ALL PLOTS	260.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	73.999231	S06N30W05	MULTIPOLYGON (((-5749.322 -376395.315, -5747.1...
55	3791777	PARKER STATION VINEYARDS-COASTAL VINEYAR	CAMP 4 - ALL PLOTS	260.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	4.678652	S06N30W05	MULTIPOLYGON (((-4980.851 -376457.889, -4979.8...
56	3790743	PARKER STATION VINEYARDS-COASTAL VINEYAR	CAMP 4 - ALL PLOTS	260.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	67.589865	S06N30W05	MULTIPOLYGON (((-4819.600 -376409.399, -4818.9...
57	3793457	PARKER STATION VINEYARDS-COASTAL VINEYAR	CAMP 4 - ALL PLOTS	260.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	94.351617	S06N30W05	MULTIPOLYGON (((-5341.376 -376455.746, -5340.7...
80	4853070	CHUMASH VINEYARDS LLC- CVCA	CAMP 4 -	280.0	ACRES	VERTEBRATE CTRL, GRAPE, WINE, UNCULTIVATED AG	277.792015	S06N30W05	MULTIPOLYGON (((-5734.564 -375663.536, -5340.6...
108	5025884	THE DAIRY -COASTAL VINEYARD CARE	THE DAIRY	116.0	ACRES	GRAPE, WINE, UNCULTIVATED AG	101.328903	S06N30W05	MULTIPOLYGON (((-6323.476 -375668.818, -6143.7...
114	5118136	THE DAIRY -CVCA	THE DAIRY-grapes	116.0	ACRES	GRAPE, WINE, UNCULTIVATED AG	101.328901	S06N30W05	MULTIPOLYGON (((-6323.476 -375668.818, -6143.7...
164	5894965	CHUMASH VINEYARDS LLC- CVCA	CAMP 4 -	280.0	ACRES	GRAPE, WINE	277.792015	S06N30W05	MULTIPOLYGON (((-5734.564 -375663.536, -5340.6...
170	5906448	RANCHO CANADA - CVCA	RANCHO CANADA-	160.0	ACRES	GRAPE, WINE, UNCULTIVATED AG	152.072178	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...
179	5954970	THE DAIRY VINEYARD-CVCA	THE DAIRY-grapes	116.0	ACRES	GRAPE, WINE, UNCULTIVATED AG	101.328901	S06N30W05	MULTIPOLYGON (((-6323.476 -375668.818, -6143.7...
227	6259346	RANCHO CANADA - CVCA	RANCHO CANADA-	160.0	ACRES	GRAPE, WINE, UNCULTIVATED AG	152.072178	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...
258	6639105	RANCHO CANADA	RANCHO CANADA-	160.0	ACRES	GRAPE, WINE, UNCULTIVATED AG	152.072178	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...
321	7295399	RANCHO CANADA	RANCHO CANADA-	160.0	ACRES	GRAPE, WINE	152.072178	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...
335	7469867	RANCHO CANADA	RANCHO CANADA-	190.0	ACRES	GRAPE, WINE	190.533840	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...
361	7848911	CHUMASH VINEYARDS LLC	CAMP 4 -	280.0	ACRES	GRAPE, WINE	277.792015	S06N30W05	MULTIPOLYGON (((-5734.564 -375663.536, -5340.6...
366	7901344	RANCHO CANADA	RANCHO CANADA-	190.0	ACRES	GRAPE, WINE	228.663901	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...
423	8534094	RANCHO CANADA	RANCHO CANADA-	190.0	ACRES	GRAPE, WINE	228.663901	S06N30W32	MULTIPOLYGON (((-6540.298 -374049.081, -6274.3...

vineyard_plot = vineyard_polygons.to_crs("EPSG:4326").hvplot.polygons(
    geo=True,
    tiles=BASEMAP,
    alpha=FIG_ALPHA,
    width=FIG_WIDTH_PX,
    height=FIG_HEIGHT_PX,
    fill_color = 'none',
    line_color = 'r',
    line_width = 3,
    title="Vineyards of Acreage greater than IQR3"
    )

hv.extension('bokeh')
vineyard_plot

Unable to display output for mime type(s): 

Interactively mapping our data

Now that we have selected a subset of the agriculture in the area we can subset our data using the rioxarray function clip.

NDVI_vineyard = NDVI.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)
NDVI_vineyard.rio.to_raster(join(data_dir,"NDVI_vineyard.tif"))
NDVI_vineyard = rxr.open_rasterio(join(data_dir,"NDVI_vineyard.tif")).squeeze("band", drop=True)
NDVI_vmin, NDVI_vmax = np.nanquantile(np.array(NDVI_vineyard), [0.02, 0.98])

ST_vineyard = ST_C.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)
ST_vineyard.rio.to_raster(join(data_dir,"ST_vineyard.tif"))
ST_vineyard = rxr.open_rasterio(join(data_dir,"ST_vineyard.tif")).squeeze("band", drop=True)
ST_vmin, ST_vmax = np.nanquantile(np.array(ST_vineyard), [0.02, 0.98])

EWT_vineyard = EWT.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)
EWT_vineyard.rio.to_raster(join(data_dir,"EWT_vineyard.tif"))
EWT_vineyard = rxr.open_rasterio(join(data_dir,"EWT_vineyard.tif")).squeeze("band", drop=True)
EWT_vmin, EWT_vmax = np.nanquantile(np.array(EWT_vineyard), [0.02, 0.98])

ET_vineyard = ET.rio.clip(vineyard_polygons.geometry.values,vineyard_polygons.crs, all_touched=True)
ET_vineyard.rio.to_raster(join(data_dir,"ET_vineyard.tif"))
ET_vineyard = rxr.open_rasterio(join(data_dir,"ET_vineyard.tif")).squeeze("band", drop=True)
ET_vmin, ET_vmax = np.nanquantile(np.array(ET_vineyard), [0.02, 0.98])

5.5 Interactive Maps

We can use holoviews to plot an interactive clickable map of these subsets over our ROI. When you click a point on the map, you can see the coordinates and data for that pixel.

clicked_values = []

# Reproject raster data
NDVI_vineyard_projected = NDVI_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)
ST_vineyard_projected = ST_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)
EWT_vineyard_projected = EWT_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)
ET_vineyard_projected = ET_vineyard.rio.reproject(WEBMAP_PROJECTION, nodata=np.nan)

# Load geospatial data
vineyard_polygons = gpd.read_file(join(data_dir,"vineyard_polygons_filtered.geojson"))

# Set up the SingleTap stream
stream = SingleTap()

# Define a function to process clicks
def interactive_click(x, y):
    if None not in [x, y]:
        NDVI_value = NDVI_vineyard_projected.sel(x=x, y=y, method="nearest").values
        ST_value = ST_vineyard_projected.sel(x=x, y=y, method="nearest").values
        ET_value = ET_vineyard_projected.sel(x=x, y=y, method="nearest").values
        EWT_value = EWT_vineyard_projected.sel(x=x, y=y, method="nearest").values
        print(x, y, NDVI_value, ST_value)
        clicked_values.append({"NDVI": NDVI_value, "ST": ST_value, "ET": ET_value, "EWT": EWT_value})

    return hv.Points((x, y)).opts(color='green', size=10)

# Create two separate plots and share the stream between them
NDVI_vineyard_map = NDVI_vineyard_projected.hvplot.image(
    # geo=True,
    rasterize=True, 
    cmap=NDVI_gradient, 
    alpha=FIG_ALPHA, 
    # tiles=BASEMAP,
    clim=(NDVI_vmin, NDVI_vmax),
    title=NDVI_title
) * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(
    line_color=LINE_COLOR,
    line_width=LINE_WIDTH,
    fill_color=FILL_COLOR
).opts(width=400, height=400)

ST_vineyard_map = ST_vineyard_projected.hvplot.image(
    # geo=True,
    rasterize=True, 
    cmap=ST_gradient, 
    alpha=FIG_ALPHA, 
    # tiles=BASEMAP,
    clim=(ST_vmin, ST_vmax),
    title=ST_title
) * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(
    line_color=LINE_COLOR,
    line_width=LINE_WIDTH,
    fill_color=FILL_COLOR
).opts(width=400, height=400)

# Duplicate the plots for the second row
EWT_vineyard_map = EWT_vineyard_projected.hvplot.image(
    # geo=True,
    rasterize=True, 
    cmap=EWT_gradient, 
    alpha=FIG_ALPHA, 
    # tiles=BASEMAP,
    clim=(EWT_vmin, EWT_vmax),
    title=EWT_title
) * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(
    line_color=LINE_COLOR,
    line_width=LINE_WIDTH,
    fill_color=FILL_COLOR
).opts(width=400, height=400)

ET_vineyard_map = ET_vineyard_projected.hvplot.image(
    # geo=True,
    rasterize=True, 
    cmap=ET_gradient, 
    alpha=FIG_ALPHA, 
    # tiles=BASEMAP,
    clim=(ET_vmin, ET_vmax),
    title=ET_title
) * vineyard_polygons.to_crs(WEBMAP_PROJECTION).hvplot.polygons(
    line_color=LINE_COLOR,
    line_width=LINE_WIDTH,
    fill_color=FILL_COLOR
).opts(width=400, height=400)

# Create a DynamicMap for interactive clicks
click = hv.DynamicMap(interactive_click, streams=[stream])

# Display the plots in a 2x2 grid with two rows along with the interactive click behavior
layout = (click * NDVI_vineyard_map * click + click * ST_vineyard_map * click +
          click * EWT_vineyard_map * click + click * ET_vineyard_map * click).cols(2)
layout

Unable to display output for mime type(s): 

Save the data from our clicked points into a csv that we can use for plotting.

# Save the clicked values to a CSV file after all clicks are processed
if clicked_values:
    df = pd.DataFrame(clicked_values)
    df.to_csv(join(data_dir,"clicked_values.csv"), index=False)
    print("Clicked values saved to 'clicked_values.csv'")
else:
    print("No clicked values to save.")

Clicked values saved to 'clicked_values.csv'

5.6 Scatter Plots

Now we can use the data saved into our csv to create scatter plots to visualize the correlations among ET, ST, EWT, and NDVI.

# Replace "your_data.csv" with the actual CSV file containing your data
df = pd.read_csv(join(data_dir,"clicked_values.csv"))

# Get combinations of all pairs of columns
combinations = [(col1, col2) for i, col1 in enumerate(df.columns) for col2 in df.columns[i + 1:]]

# Set up subplots without equal aspect ratio and sharex/sharey set to False
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(10, 12), sharex=False, sharey=False)
fig.subplots_adjust(hspace=0.5)

# Iterate over combinations and create scatter plots with best-fit lines
for (col1, col2), ax in zip(combinations, axes.flatten()):
    sns.regplot(x=col1, y=col2, data=df, ax=ax, ci=None, line_kws={"color": "red"})

    # Filter out rows with missing values in the selected columns
    filtered_df = df[[col1, col2]].dropna()

    # Calculate and plot the best-fit line using linear regression
    slope, intercept, r_value, p_value, std_err = linregress(filtered_df[col1], filtered_df[col2])
    x_range = np.linspace(filtered_df[col1].min(), filtered_df[col1].max(), 100)
    y_fit = slope * x_range + intercept
    ax.plot(x_range, y_fit, color='blue', linestyle='--', label=f'Best Fit (R={r_value:.2f})')

    ax.legend()
    ax.set_xlabel(col1)
    ax.set_ylabel(col2)
    ax.set_title(f"{col1} vs {col2}")

    # Set x and y axis ticks with different values
    ax.set_xticks(np.linspace(filtered_df[col1].min(), filtered_df[col1].max(), 5))
    ax.set_yticks(np.linspace(filtered_df[col2].min(), filtered_df[col2].max(), 5))

# Adjust layout
plt.tight_layout()

# Show the plot
plt.show()